'''
Created on Oct 7, 2010

@author: oabalbin
'''
'''
Created on Jun 15, 2010

@author: oabalbin
'''
import os
import sys
import copy
import pickle
import numpy as np
from optparse import OptionParser
from datetime import datetime
from collections import defaultdict, deque

import RNAseq.array as seqarray
import RNAseq.sampleinfo.sample_info as si
import RNAseq.arraytools.subset_genes as sg
import RNAseq.arraytools.copa as ac


if __name__ == '__main__':
    
    optionparser = OptionParser("usage: %prog [options] ")
    optionparser.add_option("-f", "--annotFile", dest="annotFile",
                            help="annotation file for all files to use")
    optionparser.add_option("-d", "--outfolder", dest="outfolder",
                            help="output folder")
    optionparser.add_option("-s", "--samplesFile", dest="samplesFile",
                            help="annotation file with the sample information")
    optionparser.add_option("-g", "--genesFolder", dest="genesFolder", default="use_all_genes",
                            help="annotation file with genes of interest")    
    optionparser.add_option("-t", "--tissueName", dest="tissueName",
                            help="name of the tissue analyzed")
    optionparser.add_option("-p", "--sampleType", dest="sampleType",
                            help="sample type. could be Cell Line or Tissue")


    (options, args) = optionparser.parse_args()
    
    
    t = datetime.now()
    tstamp = t.strftime("%Y_%m_%d_%H_%M")
    
    outfiles = ['permutated_scores']
    time_stamp_outfiles=[]
    for name in outfiles: 
        time_stamp_outfiles.append( options.outfolder+tstamp+'_'+options.sampleType+'_'+options.tissueName+'_'+name+'.txt')
        
    thisfile = open(options.annotFile)
    outfile_outliersinsample= open(time_stamp_outfiles[0],'w')
    array_dump_file = options.annotFile + '_dump_file'

    ####################
    ##### Parameters
    percentils_list =  [80,85,90,95,98] #[75,80,85,90] 
    method='mean'               # median, mean, sum
    nout2report=10
    expression_cutoff=0.0       # minimum median RPKM of expression
    base_expression_cutoff= 50  # In percentil. It is used to report the expression value of a particular gen in a sample. 
                                # The median usually can be used, however because in many cases median=1, It is maybe desireable to use for e.g the75%. 
    median_shift=1.0
    sdv_factor=0.25             #0.25
    not_report_isoforms=True    # Report gene outlier isoforms
    nrandom_permutations = 1000    # Number of random permutations
    
    
    # create the meme array
    filename=options.annotFile
    matrixfilename=options.annotFile+'memmap'
    header_cols=2
    matobj = seqarray.read_matrix(filename, matrixfilename, header_cols)

    all_samples, all_genes, all_gene_intervals = \
    seqarray.create_matobj_dicts(matobj.cols,matobj.sample_names)
    
    ####################
    # samples and gene list to use
    select_samples_by = [options.sampleType.replace('_',' '),options.tissueName,'PASS']

    sampleInfoFile = open(options.samplesFile)
    sample_info_dict = si.get_sample_info(sampleInfoFile,select_samples_by)
    sample_dict = si.generate_sample_lists(sample_info_dict)
    sample_subset, sample_cols, benign_cols = si.get_sample_columns(sample_dict, all_samples)
    tumor_cols = list(set(sample_subset.keys()).difference(benign_cols))

    if options.genesFolder!='use_all_genes':
        gene_list_file = options.genesFolder
        use_all_genes=False
    else:
        use_all_genes=True
    
    
    ########### Program               
    outfile_outliers_list = open(time_stamp_outfiles[0],'w')
    # subset of genes and samples to consider
    thisgenelist = sg.list_of_names(open(gene_list_file))
    # Get the subset of genes
    gene_subset, gene_rows, gene_subset_intervals = sg.get_gene_subset_mod(thisgenelist, all_genes, all_gene_intervals)            
    # create a sub array  using the genelist and the samplelist
    expression_matrix = matobj.exprs[gene_rows, :]
    expression_matrix = expression_matrix[:,sample_cols]
    
    thisGeneArray = sg.create_thisGenearray(expression_matrix, gene_subset, sample_subset, \
                                            list(benign_cols), gene_subset_intervals)
            
    thisGeneArray.copa_norm_across_samples(median_shift)
    outliers_cutoff = thisGeneArray.expVal.shape[0]
            
    ### Run permutation routine
    random_copa_scores = ac.copa_permute_matrix_rows(expression_matrix, benign_cols, tumor_cols, nrandom_permutations, \
                                                     percentils_list, expression_cutoff, sdv_factor, method)
    
    np.savetxt(time_stamp_outfiles[0],random_copa_scores, fmt='%10.1f', delimiter='\t')
    